% This script defines the objective function for minimum distance
% estimation of link formation probability in Poisson randome network 
% Date: April 3, 2019
function f = fitObjFun_vp(p,beta,gamma,lambda,mu,n,lb,ub,ind,l,wm);
%% summary of inputs & outputs
% INPUTS
% -- p: link formation prob in [0,1]^3 (low, medium, high, monotonic in age diff)
% -- beta, gamma, lambda: k-by-1, k-by-1 and 1-by-1 vectors of structural parameters
% -- n:  class size;  lb < n is # of rows used for fitting
% -- mu: lb-by-n-by-k array of reduced-form coefficients, where k is the #
%        of individual regressors
% -- wm: weight matrix in minimum distance estimation
% OUTPUTS
% -- f: non-negative scalar value summarizing distance between 'mu' and
%       structural estimates
%% calculate E(M) and E(M*G) under 'p' through simulation 
rng(5, 'twister'); s  = 1e4; % # of simulated networks 
rg = zeros(n,n,s); % ki = 3;
rg(ind{1},ind{1},:) = binornd(1,p(1),[l(1) l(1) s]);
rg(ind{2},ind{2},:) = binornd(1,p(1),[l(2) l(2) s]);
rg(ind{3},ind{3},:) = binornd(1,p(1),[l(3) l(3) s]);
rg(ind{1},ind{2},:) = binornd(1,p(2),[l(1) l(2) s]);
rg(ind{2},ind{1},:) = binornd(1,p(2),[l(2) l(1) s]);
rg(ind{2},ind{3},:) = binornd(1,p(2),[l(2) l(3) s]);
rg(ind{3},ind{2},:) = binornd(1,p(2),[l(3) l(2) s]);
rg(ind{1},ind{3},:) = binornd(1,p(3),[l(1) l(3) s]);
rg(ind{3},ind{1},:) = binornd(1,p(3),[l(3) l(1) s]);
rg = rg.*repmat(ones(n,n)-eye(n),[1,1,s]);
m  = zeros(n,n,s); mg = zeros(n,n,s);
for si = 1:s, % loop in simulated networks
    rgs = rg(:,:,si); 
    while min(sum(rgs,2)) == 0, 
        % draw the row again if network has "isolated" individual
        ri = find(sum(rgs,2)==0);
        for r = ri',
            rgs(r,:) = binornd(1,p(1),[1 n]);
            rgs(r,r) = 0;
        end;
    end; 
    gs  = rgs./repmat(sum(rgs,2),1,n);% row-normalization
    ms  = (eye(n) - lambda*gs)\eye(n);
    m(:,:,si)  = ms; 
    mg(:,:,si) = ms*gs;    
end;
Em  = mean(m,3);  Emg = mean(mg,3);

%% calculate weighted distance between reduced-form coefficients and implied moments
diff = [];
for ki = 1:3,
    D    =   beta(ki)*Em(1:lb,:) + gamma(ki)*Emg(1:lb,:) - mu(:,:,ki);
    diff = [ diff, mean(diag(D)) , mean(mean(D.*(ones(lb,ub)-[eye(lb) zeros(lb,ub-lb)]))) ];  
    % avg individual and contextual effects
end;
f    = diff*wm*diff'; % weighted measure of distance